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Abstract - We demonstrate that extending the Shadow Wave Function to fermionic systems 
facilitates to accurately calculate strongly-correlated multi-reference systems such as the stretched 
H 2 molecule. This development considerably extends the scope of electronic structure calculations 
and enables to efficiently recover the static correlation energy using just a single Slater determinant. 


Introduction. — One of the most outstanding prob¬ 
lems of computational physics and quantum chemistry is 
the ability to devise a quantitatively precise, yet computa¬ 
tionally tractable, method to accurately break a chemical 
bond across an entire reaction coordinate. A particularly 
simple example is the H 2 molecule, in particular when the 
covalent bond between the H atoms is stretched. Effec¬ 
tive single-particle theories, such as the widely employed 
Hartree-Fock (HF) or Density Functional Theory (DFT) 
methods, describe the covalent bond well, but the energy 
is severely overestimated upon dissociation This well- 
known problem is attributed to the multi-reference char¬ 
acter of the stretched H 2 molecule, or static electron cor¬ 
relation that arises in situations with degeneracy or near¬ 
degeneracy, as in transition metal chemistry and strongly- 
correlated systems in general [^. As a consequence, the 
stretched H 2 molecule and similar problems are typically 
dealt with using multi-determinant wave functions j^. 
However, for larger systems with many degeneracies, the 
number of determinants quickly becomes unfeasible |^[^. 

Variational Monte Carlo (VMC) is an accurate stochas¬ 
tic method that utilizes the full many-body wave func¬ 
tion (WF) and permits to approximately solve the many- 
body Schrodinger equation [^. In contrast to quantum- 
chemical methods, where the computational complexity 
grows rapidly with the number of considered electrons N 


j^, the formal scaling of quantum Monte Carlo (QMC) 
methods is similar to those of HF and DFT |8}|T^ . How¬ 
ever, since it typically relies on HF or DFT orbitals to 
construct the Slater determinant (SD), it only allows to 
extract the vast majority of dynamic electron correlation, 
but suffers from exactly the same static correlation error 
that is characteristic for single-reference electronic struc¬ 
ture methods [T^ . 

In the present work we demonstrate that extending the 
Shadow Wave Function (SWF), which was first introduced 
by Kalos and coworkers [l4)[T^ , to fermionic systems al¬ 
lows bypass the static correlation problem and permits to 
study strongly-correlated multi-reference systems within a 


much more efficient single-determinant scheme 16 -20 


Method. — Let us begin by very briefly reviewing the 
basic principles underlying the VMC method that relies 
on the Rayleigh-Ritz variational principle and importance 
sampled Monte Carlo (MC) techniques to efficiently eval¬ 
uate high-dimensional integrals [^ -23 . Due to the fact 
that the exact WF of the electronic ground state is un¬ 
known from the outset, it is approximated by a trial WF 
ipT{R,a), where R = (ri, r 2 ,... rjv) are the particle co¬ 
ordinates. The variational parameters a = 
which corresponds to the lowest variational energy 


E = 


J dRif^{R,a)IIifT{R,ot) 
f dR'tp^{R, a)'ipT{R, ct) 


( 1 ) 
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represents the best possible approximation of the elec¬ 
tronic ground state within the given trial WF 
Thus, the accuracy of a VMC simulation critically de¬ 
pends on how well the particular trial WF mimics the 
exact ground state WF. 

In order to efficiently evaluate the high-dimensional in¬ 
tegral of Eq. Q, it is conveniently rewritten as 


2. Sample from the trial WF, which corresponds to the 
initial variational parameters and estimate the 
n-dimensional vector by means of Eq. Then, 
normalize 


:= 




E = 


JdR\'iPT{R, 




( 2 ) 


which enables to compute it by means of MC methods. By 
sampling M points from the probability density function 


p{R) 




(3) 


using the M(RT)^ algorithm (also known as the Metropolis 
algorithm) , the variational energy can be estimated as 


1 ^ > 

i^l 

(4) 

where 


Eloc[R) — , , 

V't(R, a) 

(5) 

is the so-called local energy. 



Stochastic Reconfiguration. Having shown how the 
high-dimensional integral of Eq. Q can be efficiently com¬ 
puted using the M(RT)^ algorithm, it is still necessary de¬ 
termine the optimal variational parameters a, which min¬ 
imizes the variational energy. Here this is solved by devis¬ 
ing a modified version of the Stochastic Reconfiguration 
(SR) method from Sorella [^. The SR scheme prescribes 
that the variational parameters are varied according to 


3. A = 

4. = A 

5. i := i + 1 

6. Estimate by sampling from the trial WF with 

variational parameters and normalize it: 


:= 


|(5ab) I 


7. If * > 2 then A := A (l -I- 0.1 x 

8. If i > 3 then 


A := A 1^0.85-k 0.3 x 

9. If z > 5 then 

A := A I 0.75-k 0.5 x 


|q,(*) Q,(* 2) I 




(i+l-j) — ry(i-j) 






10. ab+i) = AfoW 


Sai = XY,fk{s ( 6 ) 

k=l 

where 

sik = {OkOi) — {Oi){Ok) 'I 

fk = {H){Ok)-{OkH) \ (7) 

Ok = af^ln(V'T) J 

and (•) = ('0 t| • IV't)- Once the direction in the variational 
parameters space that minimizes the variational energy 
has been computed, the step length A along this direction 
needs to be determined. Due to the fact that determin¬ 
ing the new direction 6 a is computational approximately 
equally expensive than calculating the the variational en¬ 
ergy, we found it convenient to start with a rather small 
value for A and continuously adjusting it on the fly during 
the optimization. Specifically, A is decreased whenever the 
search direction changes, which means that we are "mov¬ 
ing too fast", and increased otherwise. The modified SR 
algorithm then reads as follows: 

1. z = 0 


11. Repeat steps 5-10 until convergence. 

The most commonly employed trial WF to describe the 
electronic structure within VMC is the so-called Jastrow- 
Slater (JS) WF and consists of a single SD that is multi¬ 
plied by a simple correlation factor of the Jastrow form to 
recover most of the dynamic correlation effects : 

'f’.MR) = det{(j)air^p)) JeeiR) Jep{R,Q), ( 8 ) 

where a and /3 are the row and column indexes of the 
SDs for the spin-up and spin-down electrons, while cfa 
are the single-particle orbitals that are typically deter¬ 
mined by mean-field theories, such as HF or DFT. The 
Jastrow correlation factor J = e“ where u{rij) 

is a two-body pseudopotential, for the electron-electron 
and electron-proton interactions are denoted as Jee and 
Jep, respectively. Here we have implemented the Yukawa- 
Jastrow pseudopotential for both Jee and Jep, which is 
defined as 

— F T 

UYVK{r) = A -- (9) 

r 
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and where A and F are both variational parameters. The 
Yukawa-Jastrow pseudopotential is able to satisfy Kato’s 
cusp condition, since 

A 

UYUK(r) AF - ^r + 0(r^). (10) 

However, we have not exploited the cusp condition to fix 
one of the two parameters, but instead we have determined 
both of them by means of the modified SR algorithm. 

Shadow Wave Function. Any arbitrary trial WF can 
be systematically improved using the SWF formalism of 
Kalos and coworkers [^, which can be derived by ap¬ 
plying the imaginary-time propagator that projects 

i/'t / V'GS onto the ground state WF V'GS- In order to 
demonstrate this, let us decompose the trial WF into 

+ 00 

i’T = '^Cn4>n, ( 11 ) 

n—0 

where </>„ are the eigenfunctions of the Schrodinger equa¬ 
tion, i.e. F[(f>n = En<j>n for all n S N, with En indi¬ 
cating the associated energy eigenenvalues. Applying the 
imaginary-time propagator onto ifji we obtain 

+ 00 

e~^^tpT = (12) 

n—0 

The projector e~'^^ causes that all excited components are 
exponentially decaying so that eventually the ground 
state energy Eq is projected out, i.e. 

+ 00 

lim = lim Cne~'^^"(j)n oc (po- (13) 

r—>-oo r—>-oo 

n—0 

As a consequence, any arbitrary trial WF can be sys¬ 
tematically improved by 

e-^^'ipTiR) = {R\e-^^\p^T) (14a) 

= J dS{R\e-^^\S){S\P;T), (14b) 


where we have introduced an integral over a complete set 
of Dirac deltas \S) and omitted the inessential normaliza¬ 
tion factor. Assuming that r <C 1, we now use the Trotter 
formula to approximate 

^-riK+V) (15) 


where K represents the kinetic term of the Hamiltonian 
. Using the equality 


{x\e-^^\y) = (16) 

a 

^If some energy eigenvalues En are negative, the corresponding 
term is exponentially increasing instead of decaying. Nevertheless, 
it is always possible to add an appropriately chosen constant energy- 
shift to H, such that all excited components are again exponentially 
decaying. 


27 


where a is a normalization factor, the eventual expression 
for the improved trial WF reads as 


r.-TH 


iPt{R) = e 




(17) 

However, throughout our derivation we have assumed 
that T <gC 1. Hence, the imaginary-time propagation is 
rather short and the trial WF only slightly improved. For 
the purpose to elongate the propagation in imaginary-time 
and to solve the Schrodinger equation exactly, the de¬ 
scribed procedure needs to be applied repeatedly, which 
eventually results in a formalism similar to the path- 
integral approach 29 . Yet, there is no explicit impor¬ 
tance sampling in path-integral MC methods [^. There¬ 
fore, following our original intention to find an improved 
and computational efficient trial WF, we rather truncate 
the projection after one step and refine the obtained func¬ 
tional form at the variational level. In other words, instead 
of approaching the limit r —>■ 0, we instead substitute r by 
a variational parameter C in the gaussian term. Moreover, 
we interpret the exponential as the Jastrow corre¬ 

lation factor Jp{R) for the protons and likewise 
as the corresponding two-body correlation term JsiS) for 
the shadows. The identity {S\iPt) = Ut(< 5') implies that 
the original trial WF has to be evaluated on the shadow 
coordinates S = (s,, S 2 ,... r^r). The latter is particularly 
important for the term that dictates the symmetry of the 
SWF, which is a product of orbitals for a bosonic and a SD 
for a fermionic system, respectively. From this it follows 
that any trial WF '0 t can be systematically improved by 
shadow formalism. The resulting SWF then reads as 

V’swF(i?) = Jp(i?) MS)ijTiS). 

(18) 

From the discussion above it is apparent that the SWF 
can also be thought of as an one-step Variational Path 


Integral 31 

Although the implementation of the SWF for bosons is 
relatively straightforward, the extension to fermionic sys¬ 
tems is nontrivial due to the antisymmetry requirement of 
the WF to obey the Pauli exclusion principle. The natural 
way to devise an antisymmetric version of the SWF is to 
introduce a SD for each of the spins as a function of S, i.e. 
det((/)a(s^)) and det(0a(s^)), respectively. This results in 
the so-called Fermionic Shadow Wave Function (FSWF) 

V’FSWf(R) = JeeiR) Jep{R,Q) J dS Jse(5', i?) 

X Jsp(S',(5) det(^a(s^))det((/)a(s^)), (19) 

where Jse{S,R) is the electron-shadow and Jsp{S,Q) the 
shadow-proton Jastrow correlation factor (TblfTsy^. The 
FSWF, however, suffers from a sign problem |19|20| , which 
differs from the infamous fermion sign problem of projec¬ 
tion QMC methods such as Green’s function or diffusion 
MC 33 , but limits its applicability to relatively small 
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systems. A simple ansatz to bypass the sign problem is 
the Antisymmetric Shadow Wave Function (ASWF) 

V'aswf(.R) = Jee{R) Jep{R,Q)det{(j)a{r^p))det{(l>a{rjj)) 

X J dSe-^^^-^'>" JUS,R)JspiS,Q), ( 20 ) 

where det(0a(rp) and det(^Q(r^)) are SDs as a function 
of the electronic coordinates only j^. Even though the 
ASWF already includes many-body correlation effects of 
any order, the FSWF is superior since it accounts not only 
for symmetric, but moreover also for asymmetric, three- 


body and backflow correlation effects 34,35 


Application to the H2 molecule. — The effective¬ 
ness of the various SWFs by means of the VMC method is 
demonstrated on the H 2 molecule, whose Hamiltonian (in 
atomic units) reads as 


H 




E 


1 


ffi -1*2 1 


+ 


1 


|qi - q2| 


i=l,2 

i=1.2 


- qj 


( 21 ) 


where r represents the electronic coordinates and q the 
protonic coordinates. Since H includes the bare Coulomb 
potential, spin interactions are neglected. Due to the fact 
that the electrons of the H 2 molecule possess antiparal¬ 
lel spins, the SDs can be replaced by the orbitals them¬ 
selves. To that extend we have considered two possibili¬ 
ties. The first one is to use simple translational invariant 
plane waves (pw) orbitals, by setting 0=1 (since ki= 0 ). 
In this way, only the Jastrow correlation factor accounts 
for all the relevant physics. Alternatively, more accurate 
orbitals can be computed at the DFT level. Here we have 
employed the PWscf program of the Quantum Espresso 
package together with a pw cutoff of just 2 Ry, the bare 
Coulomb pseudopotential and the PBE exchange and cor¬ 
relation functional 37 . Since for the H 2 molecule the 


sign problem is irrelevant, it is possible to directly employ 
both, the ASWF and the more accurate FSWF. The spe¬ 
cific trial WFs we have considered in the present work are 
summarized in Tab. [TJ 

Results aud Discussiou. — In Fig. the H2 binding 
energy curves as obtained using the various trial WFs are 
shown together with the exact full configuration interac¬ 
tion (Cl) reference data for comparison jM 39 . The dif¬ 
ferences between the considered trial WFs and the exact 
full Cl reference is plotted in Fig. It is evident that 
although the DFT orbitals are generally superior to sim¬ 
ple pw orbitals, in either case the commonly employed 
JS trial WF substantially overestimates the potential en¬ 
ergy of the H 2 molecule upon dissociation. In contrast, the 
ASWF and FSWF are able to quantitatively reproduce the 
exact full Cl reference data. As can be seen in Fig. in 
particular for the FSWF, the energy difference is not only 



r [Bohr] 


Figure 1: Variational binding energy curves of the H 2 molecule, 
obtained by subtracting the variational energy with r = 
10 Bohr (dissociated energy) to the total variational energy. 
We remark that, as a consequence of such offset, a lower bind¬ 
ing energy does not imply a lower variational energy, since such 
offset is different for each trial wave function. 


very small, but more importantly approximately constant. 
This is to say that despite the strong multi-reference char¬ 
acter of the stretched H 2 molecule, the FSWF is not only 
capable to recover the dynamic but also the static corre¬ 
lation energy, using a single SD only. As such, even with 
a nearly minimal basis set, the FSWF is very competitive 
with exact or highly accurate, but computational much 
more demanding, electronic structure techniques such as 
full Cl QMC [^[^, projection QMC 33 and the Cou¬ 
pled Cluster method, which for the H 2 molecule is exact 
and equivalent to the full Cl method 1^. 


Conclusions. — To summarize, in the present work 
we have extended the SWF to fermionic systems and de¬ 
rived the underlying connection between the SWF formal¬ 
ism and projection QMC methods via the imaginary-time 
propagator. Moreover, we have demonstrated that the 
commonly used JS-DFT trial WF is able to accurately 
describe the covalent H 2 bond, but fails to recover most 
of the static correlation energy of the stretched dimer, 
which possess a sizable multi-reference character and typi¬ 
cally would require the usage of a multi-determinant WF. 
However, the ASWF and especially the FSWF permits 
to study strongly-correlated multi-reference systems and 
is able to quantitatively reproduce the exact H 2 binding 
energy curve within an very efficient single-determinant 
QMC method. 
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JS-pw 

Jee{R)Jep{R,Q) 

JS-DFT 

Jee{R) JepiR, Q) (r^)(()^^"^ (r'^) 

ASWF-pw/FSWF-pw 

Jee{R)Jep{R,Q) / 

dSJUS,R)Jsp{S,Q) 

ASWF-DFT 

Jee(i?) Jep{R, Q) (r^)(/)°^^ (r'^) / dS J^S, R) Jsp{S, Q) 

FSWF-DFT 

Jee{R)Jep{R,Q) ^ 

dS JUS, R) Jsp{S, Q) 


Table 1: Employed trial WFs for the unpolarized H 2 molecule. 



Figure 2: Difference between the binding energy curves of the 
H 2 molecules as obtained using the various trial WFs and the 
exact full Cl reference |38||39| . 

ing time through the John von Neumann Institute for 
Computing (NIC) on the CCS share of the supercom¬ 
puter JUQUEEN at the Jiilich Supercomputing Centre 
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